/*
	File:		Spectrum.cc

    Function:	Provides classes for manipulating spectra

	Author:		Andrew Willmott

	Notes:		

*/

#include <stdlib.h>
#include "Spectrum.h"
#include "ColourSystem.h"
#include "VecUtil.h"

const Int	Spectrum::kComponents = 100;	// number of basis functions
const Float	Spectrum::kStart = 350.0;		// start of spectrum
const Float	Spectrum::kEnd = 800.0;			// end of spectrum

// CIE colour matching functions, x_bar, y_bar and z_bar
// 360-830nm in steps of 5nm
// good source for this stuff is http://cvision.ucsd.edu/basicindex.htm

const Float kCIE_Start = 360.0;
const Float kCIE_End = 830.0;
const Float kCIE_Delta = 5.0;

static Float tCIE_xbar[95] =
{
#include "cie_xbar.spectrum"	
};
RegularSpectralCurve kCIE_xbar(tCIE_xbar, 360, 830, 95);

static Float tCIE_ybar[95] =
{
#include "cie_ybar.spectrum"	
};
RegularSpectralCurve kCIE_ybar(tCIE_ybar, 360, 830, 95);

static Float tCIE_zbar[95] =
{
#include "cie_zbar.spectrum"	
};
RegularSpectralCurve kCIE_zbar(tCIE_zbar, 360, 830, 95);

Colour Spectrum::ToCIE_XYZ()
{
	Int 	i;
	Colour	c = vl_0;
	
	for (i = 0; i < kComponents; i++)
	{
		ClrReal x = (i + 0.5) / kComponents;
		ClrReal lambda = mix(kStart, kEnd, x);
		
		c[0] += kCIE_xbar.Sample(lambda) * SELF[i];
		c[1] += kCIE_ybar.Sample(lambda) * SELF[i];
		c[2] += kCIE_zbar.Sample(lambda) * SELF[i];
	}

	return(c);
}

// --- Spectral curve methods -------------------------------------------------

SpectralCurve::~SpectralCurve()
{
}

SpectralCurve::operator Spectrum()
{
	Spectrum	result;
	Float		lambda, x;
	Int			i;

#ifdef USE_INT
	Float		delta, integral, lastIntegral;
	cout << "INTEGRATING" << endl;
	delta = (Spectrum::kEnd - Spectrum::kStart) / Spectrum::kComponents;
	lambda = Spectrum::kStart;
	lastIntegral = Integral(lambda);
	
	for (i = 0; i < Spectrum::kComponents; i++)
	{
		lambda += delta;
		integral = Integral(lambda);
		// normalised box functions, so we divide by the support size
		result[i] = (integral - lastIntegral) / delta;
		lastIntegral = integral;
	}
#else

	for (i = 0; i < Spectrum::kComponents; i++)
	{
		x = (i + 0.5) / Spectrum::kComponents;
		lambda = Spectrum::kStart + x * (Spectrum::kEnd - Spectrum::kStart);

		result[i] = Sample(lambda);
	}

#endif
	
	return(result);
}

Float SpectralCurve::Integral(Float lambda) const
{
	const Float	kIntegrateDelta = 1.0;
	Float		s, sum = 0.0;
	Int			n = 0;
	
	for (s = Spectrum::kStart; s < lambda; s += kIntegrateDelta)
	{
		sum += Sample(s);
		n++;
	}
	
	if (n > 0)
		sum *= (lambda - Spectrum::kStart) / n;

	return(sum);
}

Float SpectralCurve::MaxCmpt() const
{
	const Float	kMaxDelta = 1.0;
	Float		s, max = Sample(Spectrum::kStart);
	
	for (s = Spectrum::kStart + kMaxDelta; s <= Spectrum::kEnd; s += kMaxDelta)
	{
		max = Max(max, Sample(s));
	}
	
	return(max);
}

Void SpectralCurve::MakeImage(Image &image, ClrReal maxAmplitude)
// Plot spectrum into image
{
	Int		i, j;
	Float	*samples = new Float[image.Width()], ampMax = 0.0;
	Colour	*waveColour = new Colour[image.Width()];

	image.Clear(cBlack);

	for (i = 0; i < image.Width(); i++)
	{
		Float	x = (i + 0.5) / image.Width();
		Float	lambda = mix(Spectrum::kStart, Spectrum::kEnd, x);
		
		samples[i] = Sample(lambda);
		waveColour[i] = WavelengthToXYZ(lambda);
		waveColour[i] = csDisplay.XYZToGamutClip(waveColour[i]);
		ampMax = Max(samples[i], ampMax);
	}

	if (maxAmplitude < 0.0)
		maxAmplitude = ampMax;

	for (j = 0; j < image.Height(); j++)
	{ 
		Float y = maxAmplitude * (j + 0.5) / image.Height();
		
		for (i = 0; i < image.Width(); i++)
			if (samples[i] >= y)
				image.SetPixel(i, j, waveColour[i]);
	}

	delete[] samples;
	delete[] waveColour;
}

Colour SpectralCurve::ToCIE_XYZ()
{
	Float	lambda, sample;
	Colour	c = vl_0;
	
	for (lambda = kCIE_Start; lambda <= kCIE_End; lambda += kCIE_Delta)
	{
		sample = Sample(lambda);
		c[0] += kCIE_xbar.Sample(lambda) * sample;
		c[1] += kCIE_ybar.Sample(lambda) * sample;
		c[2] += kCIE_zbar.Sample(lambda) * sample;
	}

	return(c);
}

static Float *MakeAccumData(Float *data, Int divisions)
{
	Int		i;
	Float	*result, sum = 0;

	result = new Float[divisions];

	for (i = 0; i < divisions; i++)
	{
		sum += data[i];
		result[i] = sum;
	}

	return(result);
}

RegularSpectralCurve::RegularSpectralCurve(Float *sdata, Int start, Int end, Int divisions) : 
	data(sdata),
	waveDivs(divisions),
	waveStart(start),
	waveEnd(end)
{
	waveDelta = (end - start) / (divisions - 1);
	accumData = MakeAccumData(data, divisions);
}

RegularSpectralCurve::~RegularSpectralCurve()
{
	delete[] accumData;
}

Float RegularSpectralCurve::Sample(Float lambda) const
// interpolate between two closest samples
{
	Float 	result, x, dx;
	Int		bin1, bin2;

	if (lambda < waveStart || lambda > waveEnd)
		return(0.0);

	x = (lambda - waveStart) / waveDelta;
		
	bin1 = (Int) x;
	dx = x - bin1;
	if (dx > 1e-8)
		bin2 = bin1 + 1;
	else
		bin2 = bin1;

	CheckRange(bin1, 0, waveDivs, "Bad bin number");
	CheckRange(bin2, 0, waveDivs, "Bad bin number");

	result = mix(data[bin1], data[bin2], dx);

	return(result);
}


IrregularSpectralCurve::IrregularSpectralCurve(Float *amps, Float *wavs, Int divs) :
	amplitudes(amps),
	wavelengths(wavs),
	divisions(divs)
{
	accumAmplitudes = MakeAccumData(amps, divs);
}

IrregularSpectralCurve::IrregularSpectralCurve() :
	amplitudes(0),
	wavelengths(0),
	divisions(0),
	accumAmplitudes(0)
{
}

IrregularSpectralCurve::~IrregularSpectralCurve()
{
	delete[] accumAmplitudes;
}

IrregularSpectralCurve &IrregularSpectralCurve::ReadSDFile(FileName filename)
{
	String			line;
	StrConstArray	tokens;
	ifstream 		s;
	Int				i = 0;
	
	cout << "parsing " << filename.GetPath() << endl;

	s.open(filename.GetPath());

	if (!s)
	{
		cerr << "(ReadSDFile) Cannot access " << filename.GetPath() << endl;
		return(SELF);
	}

	while (s)
	{
		if (line.ReadLine(s))
		{
			Split(line, tokens);
			if (tokens.NumItems() == 0 || tokens[0][0] == '#')
				;
			else if (tokens.NumItems() == 1)
			{
				divisions = atoi(tokens[0]);
				cout << "setting size to " << divisions << endl;
				delete[] wavelengths;
				delete[] amplitudes;
				wavelengths = new Float[divisions];
				amplitudes = new Float[divisions];
				i = 0;
			}
			else if (tokens.NumItems() >= 2)
			{
				if (!wavelengths)
					_Error("(ReadSDFile) didn't find number of divisions");
				wavelengths[i] = atof(tokens[0]);
				amplitudes[i] = atof(tokens[1]);
				i++;
			}
		}
	}	

	s.close();
	accumAmplitudes = MakeAccumData(amplitudes, divisions);
	
	return(SELF);
}

Float IrregularSpectralCurve::Sample(Float lambda) const
// interpolate between two closest samples
// need binary search or something to speed this up. (cache?)
// currently we extend the endpoint samples out to infinity either
// side. Perhaps should ramp down over X nm instead?
{
	Float 	x, result;
	Int		i;

	if (divisions == 0)
		return(0.0);
	else if (divisions == 1)
		return(amplitudes[0]);

	if (lambda < wavelengths[0])
		return(amplitudes[0]);

	i = 1;
	while (i < divisions)
		if (lambda < wavelengths[i])
		{
			x = (lambda - wavelengths[i - 1])
					/ (wavelengths[i] - wavelengths[i - 1]);

			result = mix(amplitudes[i - 1], amplitudes[i], x);

			return(result);
		}
		else
			i++;

	return(amplitudes[divisions - 1]);
}


// --- Utility routines -------------------------------------------------------

Colour WavelengthToXYZ(ClrReal lambda)
{
	Colour c;

	// calculate XYZ
	c[0] = kCIE_xbar.Sample(lambda);
	c[1] = kCIE_ybar.Sample(lambda);
	c[2] = kCIE_zbar.Sample(lambda);
	
	return(c);
}

// chromaticity stuff originally from utah's sunsky code

// 300-830 10nm
static Float tS0Amplitudes[54] =
{
	0.04,6.0,29.6,55.3,57.3,
	61.8,61.5,68.8,63.4,65.8,
	94.8,104.8,105.9,96.8,113.9,
	125.6,125.5,121.3,121.3,113.5,
	113.1,110.8,106.5,108.8,105.3,
	104.4,100.0,96.0,95.1,89.1,
	90.5,90.3,88.4,84.0,85.1,
	81.9,82.6,84.9,81.3,71.9,
	74.3,76.4,63.3,71.7,77.0,
	65.2,47.7,68.6,65.0,66.0,
	61.0,53.3,58.9,61.9
};

static Float tS1Amplitudes[54] =
{
	0.02,4.5,22.4,42.0,40.6,
	41.6,38.0,42.4,38.5,35.0,
	43.4,46.3,43.9,37.1,36.7,
	35.9,32.6,27.9,24.3,20.1,
	16.2,13.2,8.6,6.1,4.2,
	1.9,0.0,-1.6,-3.5,-3.5,
	-5.8,-7.2,-8.6,-9.5,-10.9,
	-10.7,-12.0,-14.0,-13.6,-12.0,
	-13.3,-12.9,-10.6,-11.6,-12.2,
	-10.2,-7.8,-11.2,-10.4,-10.6,
	-9.7,-8.3,-9.3,-9.8
};

static Float tS2Amplitudes[54] =
{
	0.0,2.0,4.0,8.5,7.8,
	6.7,5.3,6.1,3.0,1.2,
	-1.1,-0.5,-0.7,-1.2,-2.6,
	-2.9,-2.8,-2.6,-2.6,-1.8,
	-1.5,-1.3,-1.2,-1.0,-0.5,
	-0.3,0.0,0.2,0.5,2.1,
	3.2,4.1,4.7,5.1,6.7,
	7.3,8.6,9.8,10.2,8.3,
	9.6,8.5,7.0,7.6,8.0,
	6.7,5.2,7.4,6.8,7.0,
	6.4,5.5,6.1,6.5
};

RegularSpectralCurve	kS0Spectrum(tS0Amplitudes, 300, 830, 54);
RegularSpectralCurve	kS1Spectrum(tS1Amplitudes, 300, 830, 54);
RegularSpectralCurve	kS2Spectrum(tS2Amplitudes, 300, 830, 54);

ChromaticitySpectrum::ChromaticitySpectrum(ClrReal x, ClrReal y)
{
    M1 = (-1.3515 - 1.7703 * x +  5.9114 * y) / (0.0241 + 0.2562 * x - 0.7341 * y);
    M2 = ( 0.03   -31.4424 * x + 30.0717 * y) / (0.0241 + 0.2562 * x - 0.7341 * y);
}

Float ChromaticitySpectrum::Sample(Float lambda) const
{
	return(kS0Spectrum.Sample(lambda) + M1 * kS1Spectrum.Sample(lambda)
			+ M2 * kS2Spectrum.Sample(lambda));
}

static Coord WavelengthToChroma(Float lambda)
{
	Colour c;

	// calculate XYZ
	c[0] = kCIE_xbar.Sample(lambda);
	c[1] = kCIE_ybar.Sample(lambda);
	c[2] = kCIE_zbar.Sample(lambda);

	c /= (c[0] + c[1] + c[2]);

	return(Coord(c[0], c[1]));
}

Void DrawChromaticityImage(Image &image)
// Plot CIE colour diagram
{
	Int		i, j;
	Colour	chromaClr;
	Coord	a1, a2, b1, b2;
	Float	aLambda, bLambda;
	Int		startX, endX;
	Bool	lastClip, clip;
	
	// the only tricky part to drawing the CIE
	// diagram is culling everything outside the (x,y)
	// curve defined by the pure wavelengths, and the
	// "line of purples" between the highest and lowest
	// wavelength in the spectrum.

	// to take care of this we use a1/a2/b1/b2 to step
	// through the spectral curve from either end, starting
	// from the beginning, and pick the correct start and
	// ending x values for each scanline.

	image.Clear(cGrey50);

	a1 = WavelengthToChroma(kCIE_Start);
	aLambda = kCIE_Start + kCIE_Delta;
	a2 = WavelengthToChroma(aLambda);
	b1 = a1;
	bLambda = kCIE_End;
	b2 = WavelengthToChroma(bLambda);
	
	for (j = 0; j < image.Height(); j++)
	{ 
		Float y = (j + 0.5) / image.Height();
		
		if (y < a1[1])
			continue;
		// use 'while' here because some consecutive wavelength
		// samples will have the same (x,y) values
		while (y > a2[1])
		{
			a1 = a2;
			aLambda += kCIE_Delta;
			a2 = WavelengthToChroma(aLambda);
		}
		while (y > b2[1])
		{
			b1 = b2;
			bLambda -= kCIE_Delta;
			b2 = WavelengthToChroma(bLambda);
		}
	
		if (bLambda < aLambda)
			// all done
			break;

		startX = Int(image.Width() * mix(a1[0], a2[0], (y - a1[1]) / (a2[1] - a1[1])));
		endX = Int(image.Width() * mix(b1[0], b2[0], (y - b1[1]) / (b2[1] - b1[1])));
		if (endX >= image.Width())
			endX = image.Width() - 1;

		lastClip = true;
		
		for (i = startX; i < endX; i++)
		{
			Float x = (i + 0.5) / image.Width();

			chromaClr = Colour(x, y, (1.0 - x - y));
//			chromaClr = Colour(x / y, 1.0, (1.0 - x - y) / y); 
			clip = MinCmpt(csDisplay.fromXYZ * chromaClr) < 0.0;
			if (clip != lastClip)
			{
				image.SetPixel(i, j, cYellow);
				lastClip = clip;
			}
			else
			{
				chromaClr = csDisplay.XYZToGamutClip(chromaClr);
				image.SetPixel(i, j, chromaClr);
			}
		}
	}
}

Colour CIEXYZToLuv(const Colour &cXYZ, const Colour &cWhiteXYZ)
{
	ClrReal			lumRatio = cXYZ[1] / cWhiteXYZ[1];
	const Colour	denomRat(1.0, 15.0, 3.0);
	ClrReal			denom, u, v, lumStar;
	Colour			result;
	
	if (lumRatio <= 0.008856)
		lumStar = 903.3 * lumRatio;
	else
		lumStar = 116.0 * pow(lumRatio, 1.0 / 3.0) -  16.0;

	denom = dot(denomRat, cXYZ);
	u = 4.0 * cXYZ[0] / denom;
	v = 9.0 * cXYZ[1] / denom;
	denom = dot(denomRat, cWhiteXYZ);
	u -= 4.0 * cWhiteXYZ[0] / denom;
	v -= 9.0 * cWhiteXYZ[1] / denom;

	result[0] = lumStar;
	result[1] = 13.0 * lumStar * u;
	result[2] = 13.0 * lumStar * v;

	return(result);
}

BlackBodySpectrum::BlackBodySpectrum(Float temperature)
{
	temp = temperature;
}

Float BlackBodySpectrum::Sample(Float lambda) const
{
	Float	wavelength = lambda * 1e-9;
	Float	result;

	result = 3.74183e-16 * pow(wavelength, -5.0);
	result /= exp(1.4388e-2 / (wavelength * temp)) - 1.0;

    return(result);
}
